%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fits Gaussian beam to set of beam diameters.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function BeamwaistPropagation
global  z wz lambda
figure(1);

lambda=633e-9;
z =[0 100 200 380]*1e-3;
wz=[920 980 1030 1100]*1e-6;

pa=[-100 300e-6];
pb=fminsearch(@qParameterFit,pa);

z_ =[-3000:500]*1e-3;
q0=1/(pb(1)-i*lambda/pi/pb(2)^2);
q=q0+z_;
wz_=sqrt(-lambda/pi./imag(1./q));

subplot(222),plot(z*1e2,wz*1e3,'+r',z_*1e2,wz_*1e3,'b','LineWidth',1.4)
set(gca,'xLim',[min(z_) max(z_)]*1e2,'FontSize',16,'LineWidth',1.4)
xlabel('\fontsize{16}{\itz}   (cm)');
ylabel('\fontsize{16}{\itw} ({\itz})   (mm)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Routine for das Fit-Programm. Die Funktion ist fd.
%  pa1 : Krümmungsradius
%  pa2 : Strahldurchmesser

function chi = qParameterFit(pa)
global z wz lambda

q0 =1/(pa(1)-i*lambda/pi/pa(2)^2);
q  =q0+z;
w  =sqrt(-lambda/pi./imag(1./q));
chi=sum((w-wz).^2./(1)/length(z));

subplot(222),plot(z*1e2,wz*1e3,'+r',z*1e2,w*1e3,'b')
set(gca,'FontSize',16,'LineWidth',1.4)
xlabel('\fontsize{16}{\itz}   (cm)');
ylabel('\fontsize{16}{\itw} ({\itz})   (mm)');
pause(.01)